function Reg=NewOLS(y,X);
%function [beta,u,Omega,var_beta,ttest,Rsquared]=OLS(y,X);

[T,K]=size(X);

P=inv(X'*X)*X';
z=X'*y;
P=X'*X;
beta=P\z;
%beta=P*y;%numerically the bad way of doing it....
u=y-X*beta;
%now get the variance of the residuals:
RSS=u'*u;
Omega=[u'*u]./(T-K);

%and the avriance of the OLS-estimator, var(vec(beta));

XXI=inv(X'*X);
var_beta=kron(Omega,XXI);


betadiag=zeros(size(beta));
betadiag(:)=sqrt(diag(var_beta));

ttest=beta./betadiag;


%ym=y-ones(T,1)*mean(y);
for l=1:size(y,2);
    %RSS(l)
    %RSS
    Rsquared(l,1)=1- ([RSS(l,l)]/[y(:,l)'*y(:,l)]);
    %Rsquared(l,1)= 1 - [(1-Rsquared(l,1))*((T-1)/(T-K))];
end;



Reg.b=beta;
Reg.u=u;
Reg.Omega=Omega;
Reg.tstat=ttest;
Reg.Rsquared=Rsquared;
Reg.var_b=var_beta;